#include "pathological.h"

// Ideal diffuse

void lambert_bsdf(color_t *c, material_t *m, point_t *p, vector_t *id,
	vector_t *od)
{
	if(dot(id,&p->n)*dot(od,&p->n) > 0) random_gray(c,0);
	else color_copy(c,&m->diffuse), color_scale(c,1/(2*PI));
}

float lambert_sample_bsdf(vector_t *v, material_t *m, point_t *p, vector_t *id)
{
	float r1 = ran(), r2 = ran();
	v->x = 2*cos(2*PI*r1)*sqrt(r2*(1 - r2));
	v->y = 2*sin(2*PI*r1)*sqrt(r2*(1 - r2));
	v->z = 1 - 2*r2;
	if(dot(v,&p->n)*dot(id,&p->n) > 0)
		v->x = -v->x, v->y = -v->y, v->z = -v->z;
	return 1/(2*PI);
}

// Ideal specular

void mirror_bsdf(color_t *c, material_t *m, point_t *p, vector_t *id,
	vector_t *od)
{
	vector_t v;
	mirror_sample_bsdf(&v,m,p,id);
	if(dot(&v,od) > 1 - EPSILON) color_copy(c,&m->specular);
	else random_gray(c,0);
}

float mirror_sample_bsdf(vector_t *v, material_t *m, point_t *p, vector_t *id)
{
	float vdot = -dot(&p->n,id);
	v->x = id->x + 2*vdot*p->n.x, v->y = id->y + 2*vdot*p->n.y,
		v->z = id->z + 2*vdot*p->n.z;
	normalize_vec(v);
	return 1;
}

// Ideal Fresnel reflection/refraction (dielectrics)

void fresnel_bsdf(color_t *c, material_t *m, point_t *p, vector_t *id,
	vector_t *od)
{
	int i;
	vector_t v;
	double cos1 = -dot(id,&p->n), cos2 = dot(od,&p->n),
		targetiorr = sqrt((1 - sqr(cos1))/(1 - sqr(cos2))),
		ior1, ior2;
	color_t *iorc1 = cos1 >= 0 ? refracstack->ior : &m->ior,
		*iorc2 = cos1 >= 0 || !refracstack->prev
			? &m->ior : refracstack->prev->ior;

	if(cos1*cos2 > 0) // Reflection
	{
		v.x = id->x + 2*cos1*p->n.x,
		v.y = id->y + 2*cos1*p->n.y,
		v.z = id->z + 2*cos1*p->n.z; 
		normalize_vec(&v);
		cos1 = fabs(cos1), cos2 = fabs(cos2);
		for(i = 0; i < SPEC_SAMPLES; i++)
			if(fabs(dot(&v,od)) > 1 - EPSILON)
			{
				ior1 = color_sample(iorc1,c->w[i]),
				ior2 = color_sample(iorc2,c->w[i]);
				if(ior1 > ior2 && cos1 > sqrt(1 - sqr(ior2/ior1)))
					c->i[i] = 1; // Total reflection
				else c->i[i] = (sqr((ior1*cos1 - ior2*cos2)
					/(ior1*cos1 + ior2*cos2)) + sqr((ior1*cos2 - ior2*cos1)
					/(ior1*cos2 + ior2*cos1)))/2; // Partial reflection
			}
			else c->i[i] = 0;
	}
	else // Refraction
	{
		cos1 = fabs(cos1), cos2 = fabs(cos2);
		for(i = 0; i < SPEC_SAMPLES; i++)
		{
			ior1 = color_sample(iorc1,c->w[i]),
			ior2 = color_sample(iorc2,c->w[i]);
			if(fabs(ior2/ior1 - targetiorr) < EPSILON)
				c->i[i] = 1 - (sqr((ior1*cos1 - ior2*cos2)
					/(ior1*cos1 + ior2*cos2)) + sqr((ior1*cos2 - ior2*cos1)
					/(ior1*cos2 + ior2*cos1)))/2;
			else c->i[i] = 0;
		}
	}
}

float fresnel_sample_bsdf(vector_t *v, material_t *m, point_t *p, vector_t *id)
{
	int i;
	double lambda = lambdas[(int)(ran()*(SPEC_SAMPLES - 1))],
		cos1 = -dot(id,&p->n), // Incident and reflectant angles
		abscos1 = fabs(cos1);
	color_t *iorc1 = cos1 >= 0 ? refracstack->ior : &m->ior,
		*iorc2 = cos1 >= 0 || !refracstack->prev
			? &m->ior : refracstack->prev->ior,
		iorrc;
	double ior1 = color_sample(iorc1,lambda),
		ior2 = color_sample(iorc2,lambda),
		iorr = ior1/ior2,
		cos2 = sqrt(1 - sqr(iorr)*(1 - sqr(cos1))), // Refraction angle
		preflect, ppath, piorr, x;

	if(iorr > 1 && abscos1 <= sqrt(1 - sqr(ior2/ior1))) // Total reflection
	{
		preflect = 1;
		goto reflect;
	}

	preflect = (sqr((ior1*abscos1 - ior2*cos2)/(ior1*abscos1 + ior2*cos2))
		+ sqr((ior1*cos2 - ior2*abscos1)/(ior1*cos2 + ior2*abscos1)))/2;

	if(ran() < preflect)
reflect:
		v->x = id->x + 2*cos1*p->n.x,
		v->y = id->y + 2*cos1*p->n.y,
		v->z = id->z + 2*cos1*p->n.z,
		ppath = preflect;
	else
	{
		x = (cos1 >= 0 ? 1 : -1)*(iorr*abscos1 - cos2);
		v->x = iorr*id->x + x*p->n.x,
		v->y = iorr*id->y + x*p->n.y,
		v->z = iorr*id->z + x*p->n.z;

		ppath = 1 - preflect;
	}

	color_copy(&iorrc,iorc2),
	color_rcp(&iorrc),
	color_mul(&iorrc,iorc1);
	for(i = 0, piorr = 0; i < SPEC_SAMPLES; i++)
		piorr += i == SPEC_SAMPLES - 1
			? iorr == iorrc.i[i]
				? 1 : 0
			: (iorr - iorrc.i[i])*(iorr - iorrc.i[i + 1]) <= EPSILON
			? 1.0/(1 + fabs(iorrc.i[i] - iorrc.i[i + 1])) : 0;
	piorr /= SPEC_SAMPLES;

	normalize_vec(v);

	return ppath*piorr;
}

